Let’s start by reading your data:
# Import libraries:
# Load libraries
# Note: commented-out libraries are experimental libraries for other functions not used in this final analysis file
library(limma)
#library(DESeq2)
library(edgeR)
#library(Rsubread)
library(sns)
Package: sns, Version: 1.1.2
Metropolis-Hastings MCMC using Stochastic Newton Sampler
Scientific Computing Group, Sentrana Inc. &
Imperial College London
#library(Glimma)
library(data.table)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
data.table 1.14.2 using 20 threads (see ?getDTthreads). Latest news: r-datatable.com
library(tidyr)
#library(magrittr)
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ dplyr 1.0.9
✔ tibble 3.1.7 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.2
✔ purrr 0.3.4
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::between() masks data.table::between()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::first() masks data.table::first()
✖ dplyr::lag() masks stats::lag()
✖ dplyr::last() masks data.table::last()
✖ purrr::transpose() masks data.table::transpose()
#library(sjPlot)
#library(sjlabelled)
library(sjmisc)
Learn more about sjmisc with 'browseVignettes("sjmisc")'.
Attaching package: ‘sjmisc’
The following object is masked from ‘package:purrr’:
is_empty
The following object is masked from ‘package:tibble’:
add_case
The following object is masked from ‘package:tidyr’:
replace_na
library(ggplot2)
#library(effects)
#library(broom)
#library(GGally)
#library(emmeans)
#library(ggpubr)
#theme_set(theme_pubr())
library(corncob)
library(phyloseq)
# Import data:
otus <- read.csv("gut_full.csv", check.names=F)
taxa <- read.csv("taxa.csv",sep='\t')
otus$sex <- factor(otus$sex) #factorize sex
otus$vendor_dashboard <- factor(otus$vendor_dashboard)
otus$bowel <- factor(as.numeric(factor(otus$bowel, levels = c(1,2,3,4), labels = c(1,2,3,4))))
otus <- otus[ which(otus$vendor_dashboard == "Second Genome" | otus$vendor_dashboard == "research-microbiome"),] # keep only data where vendor is explicit
otus$vendor_dashboard = str_replace_all(otus$vendor_dashboard,"research-microbiome","DNA Genotek")
otus$vendor_dashboard <- factor(otus$vendor_dashboard) # factorize vendor
#otus[ which(otus$vendor_dashboard == "Second Genome"),]
otus <- within(otus, bowel <- relevel(bowel, ref = 3))
# Pre-processing and checking validity of data
# Algorithms provided by Christian Diener, PhD:
dat <- otus
otus <- otus[!duplicated(otus$public_client_id), ]
genus_cols <- grepl("taxa_", names(otus))
rownames(otus) <- otus$public_client_id
sdata <- otus[, !genus_cols]
bowel <- paste0(sdata['bowel'])
table(sdata['bowel'])
3 1 2 4
1796 83 769 41
otus <- as.matrix(otus[, genus_cols])
colnames(otus) <- gsub("taxa_", "", colnames(otus))
tax_matrix <- as.matrix(taxa[, 2:ncol(taxa)])
rownames(tax_matrix) <- taxa[, 1]
as.data.frame(otus)
Now we convert it to fit in a phyloseq object. We need to fulfill the following rules:
To double check if everything works let’s do some validity checks:
stopifnot(all(rownames(otus) %in% rownames(sdata)))
stopifnot(nrow(otus) == nrow(sdata))
print("sample names match")
[1] "sample names match"
stopifnot(all(colnames(otus) %in% rownames(tax_matrix)))
stopifnot(ncol(otus) == nrow(tax_matrix))
stopifnot(!anyDuplicated(tax_matrix))
print("taxa look okay")
[1] "taxa look okay"
stopifnot(!anyDuplicated(sdata))
print("sample data looks okay")
[1] "sample data looks okay"
If that passes we can go ahead and build our phyloseq object.
ps <- phyloseq(
otu_table(otus, taxa_are_rows = FALSE),
tax_table(tax_matrix),
sample_data(sdata)
)
ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 68 taxa and 2689 samples ]
sample_data() Sample Data: [ 2689 samples by 7 sample variables ]
tax_table() Taxonomy Table: [ 68 taxa by 6 taxonomic ranks ]
names(sample_data(ps))
[1] "public_client_id" "sex" "age" "BMI_CALC" "vendor_dashboard" "bowel" "eGFR"
bowel <- sample_data(ps)$bowel
#Differential Test
dv_analysis <- differentialTest(formula = ~ bowel + sex + age + BMI_CALC + vendor_dashboard + eGFR,
phi.formula = ~ 1,
formula_null = ~ sex + age + BMI_CALC + vendor_dashboard + eGFR,
phi.formula_null = ~ 1,
data = ps,
test = "LRT", boot = FALSE,
full_output = TRUE,
fdr_cutoff = 0.05)
#See the significant taxa from the computation:
dv_analysis$significant_taxa
[1] "Bacteria.Actinobacteria.Actinobacteria.Bifidobacteriales.Bifidobacteriaceae.Bifidobacterium" "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Coriobacteriaceae.Collinsella"
[3] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Eggerthellaceae.Adlercreutzia" "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Bacteroidaceae.Bacteroides"
[5] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Marinifilaceae.Odoribacter" "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Rikenellaceae.Alistipes"
[7] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Tannerellaceae.Parabacteroides" "Bacteria.Firmicutes.Bacilli.Lactobacillales.Streptococcaceae.Streptococcus"
[9] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.Christensenellaceae_R-7_group" "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.nan"
[11] "Bacteria.Firmicutes.Clostridia.Clostridiales.Clostridiales_vadinBB60_group.nan" "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_AD3011_group"
[13] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_UCG-001" "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.nan"
[15] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Agathobacter" "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Blautia"
[17] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnoclostridium" "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospira"
[19] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_NK4A136_group" "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_UCG-004"
[21] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Marvinbryantia" "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.nan"
[23] "Bacteria.Firmicutes.Clostridia.Clostridiales.Peptostreptococcaceae.Romboutsia" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Anaerotruncus"
[25] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Butyricicoccus" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Candidatus_Soleaferrea"
[27] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.DTU089" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Faecalibacterium"
[29] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.GCA-900066225" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Intestinimonas"
[31] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Negativibacillus" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Oscillibacter"
[33] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_5" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_9"
[35] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_NK4A214_group" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-002"
[37] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-004" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-005"
[39] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-013" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_1"
[41] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_2" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Subdoligranulum"
[43] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.UBA1819" "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.nan"
[45] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.Holdemania" "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.nan"
[47] "Bacteria.Verrucomicrobia.Verrucomicrobiae.Verrucomicrobiales.Akkermansiaceae.Akkermansia"
saveRDS(dv_analysis, "corncob_all3.rds")
#Load the .rds file as "dv_analysis" to continue from already computed dv_analysis file from before.
#This way you don't need to run the computation all over again each time you want to run the code from here below:
library(DataCombine)
taxa<-dv_analysis$all_models
dv_analysis$significant_taxa
[1] "Bacteria.Actinobacteria.Actinobacteria.Bifidobacteriales.Bifidobacteriaceae.Bifidobacterium"
[2] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Coriobacteriaceae.Collinsella"
[3] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Eggerthellaceae.Adlercreutzia"
[4] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Bacteroidaceae.Bacteroides"
[5] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Marinifilaceae.Odoribacter"
[6] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Rikenellaceae.Alistipes"
[7] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Tannerellaceae.Parabacteroides"
[8] "Bacteria.Firmicutes.Bacilli.Lactobacillales.Streptococcaceae.Streptococcus"
[9] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.Christensenellaceae_R-7_group"
[10] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.nan"
[11] "Bacteria.Firmicutes.Clostridia.Clostridiales.Clostridiales_vadinBB60_group.nan"
[12] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_AD3011_group"
[13] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_UCG-001"
[14] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.nan"
[15] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Agathobacter"
[16] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Blautia"
[17] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnoclostridium"
[18] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospira"
[19] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_NK4A136_group"
[20] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_UCG-004"
[21] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Marvinbryantia"
[22] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.nan"
[23] "Bacteria.Firmicutes.Clostridia.Clostridiales.Peptostreptococcaceae.Romboutsia"
[24] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Anaerotruncus"
[25] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Butyricicoccus"
[26] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Candidatus_Soleaferrea"
[27] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.DTU089"
[28] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Faecalibacterium"
[29] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.GCA-900066225"
[30] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Intestinimonas"
[31] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Negativibacillus"
[32] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Oscillibacter"
[33] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_5"
[34] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_9"
[35] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_NK4A214_group"
[36] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-002"
[37] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-004"
[38] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-005"
[39] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-013"
[40] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_1"
[41] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_2"
[42] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Subdoligranulum"
[43] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.UBA1819"
[44] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.nan"
[45] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.Holdemania"
[46] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.nan"
[47] "Bacteria.Verrucomicrobia.Verrucomicrobiae.Verrucomicrobiales.Akkermansiaceae.Akkermansia"
as.data.frame(dv_analysis$p_fdr)
taxa_p <- DropNA(as.data.frame(dv_analysis$p_fdr), Var = "dv_analysis$p_fdr")
0 rows dropped from the data frame because of missing values.
taxa_p[1]
[1] 0.02120539
taxa[[1]]$coefficients
Estimate Std. Error t value Pr(>|t|)
mu.(Intercept) -3.368074680 0.182186294 -18.4869816 6.486060e-72
mu.bowel1 0.108766498 0.121122817 0.8979852 3.692741e-01
mu.bowel2 0.136368966 0.046755183 2.9166599 3.567424e-03
mu.bowel4 -0.212360465 0.178644229 -1.1887340 2.346497e-01
mu.sexW -0.212789194 0.048369795 -4.3992163 1.128862e-05
mu.age -0.016148200 0.001821779 -8.8639749 1.378229e-18
mu.BMI_CALC -0.003154924 0.003370311 -0.9360928 3.493098e-01
mu.vendor_dashboardSecond Genome 0.127694322 0.060993339 2.0935782 3.639137e-02
mu.eGFR 0.001582766 0.001243639 1.2726892 2.032388e-01
phi.(Intercept) -2.808912916 0.043294681 -64.8789385 0.000000e+00
sample_data(arivale_phylo)
Sample Data: [3653 samples by 26 sample variables]:
#Tom's code:
#Rarefy genotek dataset to even depth
#rarefy even depth is a command in phyloseq
#rarefy for alpha diversity, but not for corncob
rarefied_genotek=rarefy_even_depth(arivale_phylo, sample.size = min(sample_sums(arivale_phylo)),
rngseed = 111, replace = FALSE, trimOTUs = TRUE, verbose = TRUE)
`set.seed(111)` was used to initialize repeatable random subsampling.
Please record this for your records so others can reproduce.
Try `set.seed(111); .Random.seed` for the full vector
...
25984OTUs were removed because they are no longer
present in any sample after random subsampling
...
richness <- estimate_richness(rarefied_genotek, measures=c("Shannon","Observed"))
#Refer to Christian (phytree?) on PD whole tree (phyloseq additional PD whole tree package, check it) and add it to richness column on df.
richness$public_client_id <-sample_data(rarefied_genotek)$public_client_id
#Refer to Christian (phytree?) on PD whole tree (phyloseq additional PD whole tree package, check it) and add it to richness column on df.
library(btools)
richness$Pielou <- richness$Shannon/richness$Observed
richness
hist(richness$Pielou,
main = "Pielou's Evenness (Shannon/Observed ASV)",
xlab = "Evenness",
xlim = c(0,0.05),
breaks = 200)
hist(richness$Shannon,
main = "Shannon Diversity",
xlab = "Diversity Index",
xlim = c(0,6),
breaks = 20)
hist(richness$Observed,
main = "Observed ASVs",
xlab = "ASVs",
xlim = c(0,1000),
breaks = 30)
library(tidyverse)
library(ggpmisc)
library(quantreg)
library(broom)
library(ggbeeswarm)
library(ggsignif)
library(gginnards)
library(ggpubr)
library(ggpmisc)
library(broom.mixed)
library(scales)
#library(nlme)
sam.n <- function(x){
return(c(y = mean(x), label = length(x)))
}
combinations <- list(c("Constipation","Low Normal"),
c("Constipation","High Normal"),
c("Constipation","Diarrhea"),
c("Low Normal","High Normal"),
c("Low Normal","Diarrhea"),
c("High Normal","Diarrhea"))
dat$bowel <- factor(dat$bowel, levels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))
library(compositions)
df <- read.csv("gut_full.csv", check.names=F)
colnames(df) = gsub("nan", "Unclassified", colnames(df))
comparisons = list(c("Constipation","High Normal"),c("Low Normal","High Normal"),c("Diarrhea","High Normal"))
df$bowel <- factor(df$bowel, levels = c(1,2,3,4), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))
df_otus <- dplyr::select(df, -c("public_client_id","bowel","vendor_dashboard","sex","age","BMI_CALC","eGFR"))
df_otus <- as.data.frame(clr(as.matrix(df_otus)))
df_select <- dplyr::select(df, c(1:7))
df_otus <- cbind(df_select,df_otus)
df <- df_otus
dat <- df_otus
library(ggrepel)
sam.n <- function(x){return(c(y = mean(x), label = length(x)))}
# ggplot2 plots
library(gplots)
library(ggrepel)
library(ggbreak)
library(sommer)
#BMF df
const <- c()
const_p <- c()
low <- c()
low_p <- c()
diarrhea <- c()
diarrhea_p <- c()
taxa_names <- c()
family_names <- c()
genus_names <- c()
likelihood <- c()
adj_p <- c()
names(dv_analysis$p) = gsub("nan", "Unclassified", names(dv_analysis$p))
#dv_analysis$all_models[order(dv_analysis$p[order(names(dv_analysis$p))])])
for (i in 1:length(dv_analysis$all_models[names(dv_analysis$p)])) {
const_p[i] <- dv_analysis$all_models[[i]]$coefficients[2,4]
const[i] <- dv_analysis$all_models[[i]]$coefficients[2,1]
low_p[i] <- dv_analysis$all_models[[i]]$coefficients[3,4]
low[i] <- dv_analysis$all_models[[i]]$coefficients[3,1]
diarrhea_p[i] <- dv_analysis$all_models[[i]]$coefficients[4,4]
diarrhea[i] <- dv_analysis$all_models[[i]]$coefficients[4,1]
taxa_names[i] <- names(dv_analysis$p)[i]
family_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][5]
genus_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][6]
likelihood[i] <- dv_analysis$all_models[[i]]$logL
adj_p[i] <- dv_analysis$p_fdr[[i]]
}
genus_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][6]
p_df <- bind_cols(taxa_names,family_names,genus_names, likelihood, adj_p, const,const_p,low,low_p,diarrhea,diarrhea_p)
New names:
• `` -> `...1`
• `` -> `...2`
• `` -> `...3`
• `` -> `...4`
• `` -> `...5`
• `` -> `...6`
• `` -> `...7`
• `` -> `...8`
• `` -> `...9`
• `` -> `...10`
• `` -> `...11`
names(p_df) <- c("Genera","Family","Genus","LogL","Adj.P","Const.Beta","Const.P.Val","Low.Beta","Low.P.Val","Diarrhea.Beta","Diarrhea.P.Val")
p_df$Combined <- paste(p_df$Family,p_df$Genus)
ab <- colSums(df_otus[names(df_otus) %in% paste("taxa_",names(dv_analysis$p),sep="")])/colSums(!!df_otus[names(df_otus) %in% paste("taxa_",names(dv_analysis$p),sep="")])
p_df$Mean.Abundance <- ab[paste("taxa_",p_df$Genera,sep="")]
p_df$Const.Adj.P.Val <- p.adjust(const_p, method = "fdr", n = length(const_p))
p_df$Low.Adj.P.Val <- p.adjust(low_p, method = "fdr", n = length(low_p))
p_df$Diarrhea.Adj.P.Val <- p.adjust(diarrhea_p, method = "fdr", n = length(diarrhea_p))
p_df <- p_df[order(-p_df$Mean.Abundance),]
set <- subset(p_df[which(p_df$Adj.P < 0.05),], Genus != 'Unclassified')
set <- subset(set[order(-set$Mean.Abundance),], Genus != 'Unclassified')
abundant_list <- list()
abundant_genus <- list()
abundant <- paste("taxa_",set[order(-set$Mean.Abundance),]$Genera,sep="")
for (i in 1:68) {
abundant_list[i] <- paste(strsplit(abundant[i],'.',fixed=TRUE)[[1]][5],strsplit(abundant[i],'.',fixed=TRUE)[[1]][6])
abundant_genus[i] <- paste(strsplit(abundant[i],'.',fixed=TRUE)[[1]][6])
}
#top 10 most abundant non-unclassified hits
abundant_trunc <- set[match(abundant_list[!is.na(set$Combined[match(abundant_list,set$Combined)])],set$Combined),][which(set$Adj.P<0.05),]$Combined[1:10]
`%!in%` <- Negate(`%in%`)
t <- set[set$Combined %!in% abundant_trunc,]
#t[!is.na(t[order(t$Adj.P),]$Mean.Abundance),][order(t$Adj.P),]
span <- rbind(set[order(set[-order(set[!set$Combined %in% abundant_trunc,]$Adj.P),][which(set$Adj.P<0.05),]$Mean.Abundance %!in% abundant_trunc),][1:10,],
t[match('Akkermansiaceae Akkermansia',t$Combined),],
t[order(t$Const.Adj.P.Val),][1:9,])
span$Letters <- LETTERS[1:20]
p_df$Letters <- c(NA)
test <- rbind(span,p_df[p_df$Combined %!in% span$Combined,])
test <- test[!duplicated(test$Genera),]
#test[order(test[!duplicated(test),]$Letters),]
p_df <- test
p_df <- p_df[order(p_df$Adj.P),]
p_df <- p_df[order(-p_df$Mean.Abundance),]
model = y ~ x
df_num <- merge(richness, df_otus, by="public_client_id")
df_num <- within(df_num, bowel <- relevel(bowel, ref = 'High Normal'))
d1_plot <- lm(Pielou ~ bowel, data = df_num)
pval1 = summary(d1_plot)$coefficients[2,4]
r1 = summary(d1_plot)$adj.r.squared
pval1
[1] 0.001840045
r1
[1] 0.01204487
d1 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Pielou, group = factor(bowel))) +
scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r1),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval1),3)))))+
geom_beeswarm(aes(color = bowel), cex = 0.5) +
geom_boxplot(alpha=0) +
ylim(0.005,0.05)+
geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
ggtitle(label = "C)\nPielou's Evenness vs BMF") +
xlab("Bowel Movement Frequency (BMF)") +
ylab("Pielou's Evenness") +
scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
guides(colour = guide_legend(override.aes = list(size=7))) +
theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d1
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (position_beeswarm).
model = y ~ x
d2_plot <- lm(Shannon ~ bowel, data = df_num)
pval2 = summary(d2_plot)$coefficients[2,4]
r2 = summary(d2_plot)$adj.r.squared
pval2
[1] 0.00964909
r2
[1] 0.01384
d2 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Shannon, group = bowel)) +
scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r2),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval2),3)))))+
geom_beeswarm(aes(color = factor(bowel)), cex = 0.5) +
geom_boxplot(alpha=0) +
geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
ggtitle(label = "B)\nShannon Diversity vs BMF") +
xlab("Bowel Movement Frequency (BMF)") +
ylab("Shannon Diversity") +
scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
guides(colour = guide_legend(override.aes = list(size=7))) +
theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d2
model = y ~ x
d3_plot <- lm(Observed ~ bowel, data = df_num)
pval3 = summary(d3_plot)$coefficients[2,4]
r3 = summary(d3_plot)$adj.r.squared
pval3
[1] 5.448268e-05
r3
[1] 0.0208449
d3 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Observed, group = bowel))+
scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r3),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval3),3)))))+
geom_beeswarm(aes(color = factor(bowel)), cex = 0.5) +
geom_boxplot(alpha=0) +
geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
ggtitle(label = "A)\nObserved ASVs vs BMF") +
xlab("Bowel Movement Frequency (BMF)") +
ylab("Observed ASVs") +
scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
guides(colour = guide_legend(override.aes = list(size=7))) +
theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d3
PSO <- ggarrange(d3, d2, d1, legend = "top", align = "hv", common.legend = TRUE, widths = c(7,7,7), heights = c(20,20,20), nrow = 1, ncol = 3)+ theme(plot.margin = margin(0,0,0,0))
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (position_beeswarm).
PSO
ggsave(
"PSOvBMF.png",
plot = PSO,
device = NULL,
path = NULL,
scale = 1,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 7.29 x 4.51 in image
# Significant Genera Plot (plot 0)
plot0s <- ggplot(p_df, mapping = aes(x=Const.Beta,y = -log10(Const.P.Val), color = ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"),group = factor(ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"))))+
geom_vline(xintercept = 0)+
geom_point(size=1, aes(x = Const.Beta, y = -log10(Const.P.Val), group = factor(ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Letters,""))),position = position_dodge(width = 1)) +
geom_label_repel(aes(x = Const.Beta, y = -log10(Const.P.Val), group = factor(ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Letters,""))), label.size = 0.1, label.padding = 0.1, box.padding = 0.5, min.segment.length = 0, point.size = NA, position = position_dodge(width = 1), size=2,label = factor(ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Letters,"")),color = "black", check_overlap = FALSE, max.overlaps = 127)+
#ggtitle("Significant Genera") +
xlab(bquote(beta~" Coefficient")) +
ylab(bquote("-log"[10]~"(P value)")) +
#scale_x_break(c(-0.6, -18)) +
scale_x_continuous(name = bquote(atop(beta["BMF"]~" Coefficient",italic("Constipation"))),
guide = guide_axis(n.dodge = 2), limits = c(-1.5,1.5)) +
theme(text =
element_text(size = 14),
plot.title = element_text(vjust = 0.5),
#plot.subtitle = element_text(size=8, hjust = 0.5),
legend.title = element_blank(),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 8),
legend.text = element_text(size = 3),
legend.position = "top",
legend.key.size = unit(0.1, "lines"))+
guides(shape = guide_legend(override.aes = list(size = 1.5)),
color = guide_legend(override.aes = list(size = 1.5), nrow = 25),
fill=guide_legend(title=NULL),
aspect.ratio = 0.95)
Warning: Ignoring unknown parameters: check_overlap
plot0s
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
Warning: Removed 13 rows containing missing values (geom_label_repel).
# Significant Genera Plot (plot 0)
plot0m <- ggplot(p_df, mapping = aes(x=Low.Beta,y = -log10(Low.P.Val), color = ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"),group = factor(ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"))))+
geom_vline(xintercept = 0)+
geom_point(size=1, aes(x = Low.Beta, y = -log10(Low.P.Val), group = factor(ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Letters,""))),position = position_dodge(width = 1)) +
geom_label_repel(aes(x = Low.Beta, y = -log10(Low.P.Val), group = factor(ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Letters,""))), label.size = 0.1, label.padding = 0.1, box.padding = 0.5, min.segment.length = 0, point.size = NA, position = position_dodge(width = 1), size=2,label = factor(ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Letters,"")),color = "black", check_overlap = FALSE, max.overlaps = 127)+
#ggtitle("Significant Genera") +
xlab(bquote(beta~" Coefficient")) +
ylab(bquote("-log"[10]~"(P value)")) +
#scale_x_break(c(-0.6, -18)) +
scale_x_continuous(name = bquote(atop(beta["BMF"]~" Coefficient",italic("Low Normal"))),
guide = guide_axis(n.dodge = 2)) +
theme(text =
element_text(size = 14),
plot.title = element_text(vjust = 0.5),
#plot.subtitle = element_text(size=8, hjust = 0.5),
legend.title = element_blank(),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 8),
legend.text = element_text(size = 3),
legend.position = "top",
legend.key.size = unit(0.1, "lines"))+
guides(shape = guide_legend(override.aes = list(size = 1.5)),
color = guide_legend(override.aes = list(size = 1.5), nrow = 25),
fill=guide_legend(title=NULL),
aspect.ratio = 0.95)
Warning: Ignoring unknown parameters: check_overlap
plot0m
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
Warning: Removed 27 rows containing missing values (geom_label_repel).
# Significant Genera Plot (plot 0)
plot0d <- ggplot(p_df, mapping = aes(x=Diarrhea.Beta,y = -log10(Diarrhea.P.Val), color = ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"),group = factor(ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"))))+
geom_vline(xintercept = 0)+
geom_point(size=1, aes(x = Diarrhea.Beta, y = -log10(Diarrhea.P.Val), group = factor(ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Letters,""))),position = position_dodge(width = 1)) +
geom_label_repel(aes(x = Diarrhea.Beta, y = -log10(Diarrhea.P.Val), group = factor(ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Letters,""))), label.size = 0.1, label.padding = 0.1, box.padding = 0.5, min.segment.length = 0, point.size = NA, position = position_dodge(width = 1), size=2,label = factor(ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Letters,"")),color = "black", check_overlap = FALSE, max.overlaps = 127)+
#ggtitle("Significant Genera") +
xlab(bquote(beta~" Coefficient")) +
ylab(bquote("-log"[10]~"(P value)")) +
#scale_x_break(c(-0.6, -18)) +
scale_x_continuous(name = bquote(atop(beta["BMF"]~" Coefficient",italic("Diarrhea"))),
guide = guide_axis(n.dodge = 2)) +
theme(text =
element_text(size = 14),
plot.title = element_text(vjust = 0.5),
#plot.subtitle = element_text(size=8, hjust = 0.5),
legend.title = element_blank(),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 8),
legend.text = element_text(size = 3),
legend.position = "top",
legend.key.size = unit(0.1, "lines"))+
guides(shape = guide_legend(override.aes = list(size = 1.5)),
color = guide_legend(override.aes = list(size = 1.5), nrow = 25),
fill=guide_legend(title=NULL),
aspect.ratio = 0.95)
Warning: Ignoring unknown parameters: check_overlap
plot0d
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
library(ggpubr)
figure1 <- ggarrange(plot0s,plot0m,plot0d, common.legend = FALSE, nrows = 1, widths=1,heights=3,ncols = 3, align = "hv", legend = "top") + theme(plot.margin = margin(0,0,0,0))
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
Warning: Removed 13 rows containing missing values (geom_label_repel).
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
Warning: Removed 27 rows containing missing values (geom_label_repel).
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
ggsave(
"BMFvsGenera1.png",
plot = figure1,
device = NULL,
path = NULL,
scale = 1.75,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 12.8 x 7.88 in image
figure1
Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto), :
Viewport has zero dimension(s)
sig = function(x){
if(x < 0.001){"***"}
else if(x < 0.01){"**"}
else if(x < 0.05){"*"}
else{NA}}
df <- read.csv("gut_full.csv", check.names=F)
colnames(df) = gsub("nan", "Unclassified", colnames(df))
comparisons = list(c("Constipation","High Normal"),c("Low Normal","High Normal"),c("Diarrhea","High Normal"))
df$bowel <- factor(df$bowel, levels = c(1,2,3,4), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))
df_otus <- dplyr::select(df, -c("public_client_id","bowel","vendor_dashboard","sex","age","BMI_CALC","eGFR"))
df_otus <- as.data.frame(clr(as.matrix(df_otus)))
df_select <- dplyr::select(df, c(1:7))
df_otus <- cbind(df_select,df_otus)
df <- df_otus
dat <- df_otus
#Gut Genera
df$p.value <- NA
df$pval <- NULL
df$p.val <- NULL
counter <<- 1
test = function(x,y,z,j) {
work <- NULL
results <- NULL
x = NULL
y = NULL
if (counter == 1) {z = comparisons[[1]][1]}
else if (counter == 2) {z = comparisons[[2]][1]}
else if (counter > 2) {z = comparisons[[3]][1]}
temp <- data.frame("name" = NA, "p.value" = NA, "bowel" = NA)
for (i in 1:20) {
temp$name <- paste("genus_",span$Genera[[i]],sep="")
temp$p.value <- span$Const.Adj.P.Val[[i]]
temp$bowel <- "Constipation"
work <- rbind(work,temp)
temp$p.value <- span$Low.Adj.P.Val[[i]]
temp$bowel <- "Low Normal"
work <- rbind(work,temp)
temp$p.value <- span$Diarrhea.Adj.P.Val[[i]]
temp$bowel <- "Diarrhea"
work <- rbind(work,temp)
}
if (z[[1]][1] == comparisons[[1]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1])}
else if (z[[1]][1] == comparisons[[2]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1+1]) }
else if (z[[1]][1] == comparisons[[3]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1+2])}
counter<<-counter+1
if (counter > 3) {counter<<-1}
return(results)
}
myplots <- list() # new empty list
df_test <- df
df_test[df_test==0] <- NA
for (genera in 1:20) {
y_name <- paste("taxa_",span$Genera[genera],sep="")
myplots[[genera]] <- local({
plotlim_lower = min(
df_test[!is.na(df_test[y_name]),][y_name])
plotlim_upper = max(
df_test[!is.na(df_test[y_name]),][y_name])
plotlim_bar = plotlim_lower - 4
plotlim_margin = 6
genera <- genera
plt <- ggplot(data = df_test, aes(x = bowel, y = .data[[y_name]], group = bowel)) +
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
geom_jitter(aes(color = bowel), size = 0.1, cex = 0.05) +
geom_boxplot(data = df_test, alpha=0.0,outlier.shape = NA) +
theme(text = element_text(size = 9)) +
ggtitle(label = paste(span$Family[genera],"\n",span$Genus[genera],sep="")) +
geom_signif(data = df_test, comparisons = comparisons, map_signif_level = sig, test = 'test', test.args = list(z = comparisons, j = genera), y_position = plotlim_bar,step_increase = 0.15, size = 0.5,
textsize = 1.5,
tip_length = c(0,0)) +
coord_cartesian(ylim=c(plotlim_lower,plotlim_upper),clip="off")+
labs(color = "BMF Category", y = ifelse((genera == 1 | genera == 6 | genera == 11 | genera == 16),"CLR Abundance","")) +
guides(colour = guide_legend(override.aes = list(size=7), title.position = 'left', nrow = 1, ncol = 4)) +
theme(plot.margin = unit(c(0,0,plotlim_margin,0), "cm"),
plot.title = element_text(size=5.75),
legend.title = element_text(size=10),
plot.subtitle = element_text(size=10),
legend.text = element_text(size=7),
axis.text.x = element_blank(),
axis.text.y = element_text(size=7),
axis.title.y = element_text(size=7),
axis.title.x = element_blank(),
aspect.ratio = 0.95)
})
}
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
library(gplots)
library(ggrepel)
library(ggbreak)
#Truncate the top 20 including verukkamansia
title <- list()
for(i in seq(1:20)) { title[i] <- myplots[[c(i)]]$labels$title }
plots1 <- myplots[title %in% paste(
span$Family,
"\n",
span$Genus,sep="")][1:20]
plots_final <- plots1[!sapply(plots1,is.null)]
final1 <- ggarrange(plotlist = plots_final[1:5], labels = LETTERS[1:5], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 56 rows containing non-finite values (stat_boxplot).
Warning: Removed 56 rows containing non-finite values (stat_signif).
Warning: Removed 56 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing non-finite values (stat_signif).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 85 rows containing non-finite values (stat_boxplot).
Warning: Removed 85 rows containing non-finite values (stat_signif).
Warning: Removed 85 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
final1
final2 <- ggarrange(plotlist = plots_final[6:10], labels = LETTERS[6:10], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))
Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 98 rows containing non-finite values (stat_boxplot).
Warning: Removed 98 rows containing non-finite values (stat_signif).
Warning: Removed 98 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 607 rows containing non-finite values (stat_boxplot).
Warning: Removed 607 rows containing non-finite values (stat_signif).
Warning: Removed 607 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 4 rows containing non-finite values (stat_boxplot).
Warning: Removed 4 rows containing non-finite values (stat_signif).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 7 rows containing missing values (geom_signif).
final2
final3 <- ggarrange(plotlist = plots_final[11:15], labels = LETTERS[11:15], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))
Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 271 rows containing non-finite values (stat_boxplot).
Warning: Removed 271 rows containing non-finite values (stat_signif).
Warning: Removed 271 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 591 rows containing non-finite values (stat_boxplot).
Warning: Removed 591 rows containing non-finite values (stat_signif).
Warning: Removed 591 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 147 rows containing non-finite values (stat_boxplot).
Warning: Removed 147 rows containing non-finite values (stat_signif).
Warning: Removed 147 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 240 rows containing non-finite values (stat_boxplot).
Warning: Removed 240 rows containing non-finite values (stat_signif).
Warning: Removed 240 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
final3
final4 <- ggarrange(plotlist = plots_final[16:20], labels = LETTERS[16:20], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))
Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 40 rows containing non-finite values (stat_boxplot).
Warning: Removed 40 rows containing non-finite values (stat_signif).
Warning: Removed 40 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 331 rows containing non-finite values (stat_boxplot).
Warning: Removed 331 rows containing non-finite values (stat_signif).
Warning: Removed 331 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 680 rows containing non-finite values (stat_boxplot).
Warning: Removed 680 rows containing non-finite values (stat_signif).
Warning: Removed 680 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 20 rows containing non-finite values (stat_boxplot).
Warning: Removed 20 rows containing non-finite values (stat_signif).
Warning: Removed 20 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
final4
counter<<-1
figure_main <- ggarrange(plotlist = plots_final[1:20], labels = LETTERS[1:20], legend = "top", align = "hv", common.legend = TRUE, nrow = 4, ncol = 5)
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 56 rows containing non-finite values (stat_boxplot).
Warning: Removed 56 rows containing non-finite values (stat_signif).
Warning: Removed 56 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing non-finite values (stat_signif).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 85 rows containing non-finite values (stat_boxplot).
Warning: Removed 85 rows containing non-finite values (stat_signif).
Warning: Removed 85 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 98 rows containing non-finite values (stat_boxplot).
Warning: Removed 98 rows containing non-finite values (stat_signif).
Warning: Removed 98 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 607 rows containing non-finite values (stat_boxplot).
Warning: Removed 607 rows containing non-finite values (stat_signif).
Warning: Removed 607 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 4 rows containing non-finite values (stat_boxplot).
Warning: Removed 4 rows containing non-finite values (stat_signif).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 7 rows containing missing values (geom_signif).
Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 271 rows containing non-finite values (stat_boxplot).
Warning: Removed 271 rows containing non-finite values (stat_signif).
Warning: Removed 271 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 591 rows containing non-finite values (stat_boxplot).
Warning: Removed 591 rows containing non-finite values (stat_signif).
Warning: Removed 591 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 147 rows containing non-finite values (stat_boxplot).
Warning: Removed 147 rows containing non-finite values (stat_signif).
Warning: Removed 147 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 240 rows containing non-finite values (stat_boxplot).
Warning: Removed 240 rows containing non-finite values (stat_signif).
Warning: Removed 240 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 40 rows containing non-finite values (stat_boxplot).
Warning: Removed 40 rows containing non-finite values (stat_signif).
Warning: Removed 40 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 331 rows containing non-finite values (stat_boxplot).
Warning: Removed 331 rows containing non-finite values (stat_signif).
Warning: Removed 331 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 680 rows containing non-finite values (stat_boxplot).
Warning: Removed 680 rows containing non-finite values (stat_signif).
Warning: Removed 680 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 20 rows containing non-finite values (stat_boxplot).
Warning: Removed 20 rows containing non-finite values (stat_signif).
Warning: Removed 20 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
ggsave(
"BMFvsGenera2.png",
plot = final1,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.9 x 6.76 in image
ggsave(
"BMFvsGenera3.png",
plot = final2,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.9 x 6.76 in image
ggsave(
"BMFvsGenera4.png",
plot = final3,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.9 x 6.76 in image
ggsave(
"BMFvsGenera5.png",
plot = final4,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.9 x 6.76 in image
ggsave(
"BMFvsGeneraMain.png",
plot = final_main,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.9 x 6.76 in image
ggsave(
"BMFvsGeneraMain.png",
plot = final_main,
device = NULL,
path = NULL,
scale = 0.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 7 x 7 in image